/*------------------------------------------------------------------------
 *  Copyright (C) 2010  Luis M. de la Cruz
 *
 *  This file is part of TUNA
 *
 *  TUNA is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  TUNA is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ------------------------------------------------------------------------*/

namespace Tuna {

/*!
 ***************************************************************************
 * Interpolate an scalarfield from nodes to center of volumes (1D,2D or 3D).
 * \verbatim
       out field (volumes)               in field (nodes)

 5  x--x--+--x--+--x--+--x--x        o-----o-----o-----o-----o  4   
    |     |     |     |     |        |     |     |     |     |
 4  x  x  |  x  |  x  |  x  x        |     |     |     |     |
    |     |     |     |     |        |     |     |     |     |
    +-----+-----+-----+-----+        o-----o-----o-----o-----o  3   
    |     |     |     |     |        |     |     |     |     |
 3  x  x  |  x  |  x  |  x  x        |     |     |     |     |
    |     |     |     |     |        |     |     |     |     |
    +-----+-----+-----+-----+  <===  o-----o-----o-----o-----o  2    
    |     |     |     |     |        |     |     |     |     |
 2  x  x  |  x  |  x  |  x  x        |     |     |     |     |
    |     |     |     |     |        |     |     |     |     |
    +-----+-----+-----+-----+        o-----o-----o-----o-----o  1    
    |     |     |     |     |        |     |     |     |     |
 1  x  x  |  x  |  x  |  x  x        |     |     |     |     |
    |     |     |     |     |        |     |     |     |     |
 0  x--x--+--x--+--x--+--x--x        o-----o-----o-----o-----o  0     

    0  1     2     3     4  5        0     1     2     3     4 

 output size : (nodes_x, nodes_y)    input size : (vols_x, vols_y)

 Note: vols_? = nodes_? + 1, for ? = x, y, z \endverbatim
 ***************************************************************************
 *  \author Luis M. de la Cruz [ Sun Mar 23 15:33:45 GMT 2008 ]
 ***************************************************************************
 */
template<class Tfield>
void NumUtils::interpolateToVolumes(Tfield &out, Tfield &in) 
{
    int bi = out.lbound(firstDim) + 1, ei = out.ubound(firstDim) - 1, 
	bj = out.lbound(secondDim)+ 1, ej = out.ubound(secondDim)- 1,
	bk = out.lbound(thirdDim) + 1, ek = out.ubound(thirdDim) - 1;
    
    Range I(bi,ei), J(bj,ej), K(bk,ek);

    if (out.rank() == 1) 
	   out(I) = ( in(I-1) + in(I) ) * 0.5;
	   out(bi-1) = in(bi-1);
	   out(ei+1) = in(ei);
    if (out.rank() == 2) {
//
// Internal points    
//
	out(I,J) = ( in(I,J) + in(I-1,J) + in(I,J-1) + in(I-1,J-1) ) * 0.25;
//	
// Boundaries
//
	out(I   ,0   ) = ( in(I , 0 ) + in(I-1, 0  ) ) * 0.5;
	out(I   ,ej+1) = ( in(I , ej) + in(I-1, ej ) ) * 0.5;
	out(0   ,J   ) = ( in(0 , J ) + in(0  , J-1) ) * 0.5;
	out(ei+1,J   ) = ( in(ei, J ) + in(ei , J-1) ) * 0.5;
//	
// Corners
//
	out(0   ,0   ) = in(0 , 0 );
	out(ei+1,0   ) = in(ei, 0 );
	out(0   ,ej+1) = in(0 , ej);
	out(ei+1,ej+1) = in(ei, ej);
    }
    if (out.rank() == 3) {
//
// Internal points    
//
	out(I,J,K) = ( in(I,J  ,K  ) + in(I-1,J  ,K  ) + 
		       in(I,J-1,K  ) + in(I-1,J-1,K  ) +
		       in(I,J  ,K-1) + in(I-1,J  ,K-1) + 
		       in(I,J-1,K-1) + in(I-1,J-1,K-1) ) * 0.125;
///	
/// Boundaries
///
/// Bot
	out(I,0   ,K) = ( in(I,0,K  ) + in(I-1,0,K  ) +
			  in(I,0,K-1) + in(I-1,0,K-1) ) * 0.25;
/// Top
	out(I,ej+1,K) = ( in(I,ej,K  ) + in(I-1,ej,K  ) +
			  in(I,ej,K-1) + in(I-1,ej,K-1) ) * 0.25;
/// Left
	out(0   ,J,K) = ( in(0,J,K  ) + in(0,J-1,K  ) +
			  in(0,J,K-1) + in(0,J-1,K-1) ) * 0.25;
/// Right
	out(ei+1,J,K) = ( in(ei,J,K  ) + in(ei,J-1,K  ) +
			  in(ei,J,K-1) + in(ei,J-1,K-1) ) * 0.25;
/// Back
	out(I,J,0   ) = ( in(I,J  ,0) + in(I-1,J  ,0) +
			  in(I,J-1,0) + in(I-1,J-1,0) ) * 0.25;
/// Front
	out(I,J,ek+1) = ( in(I,J  ,ek) + in(I-1,J  ,ek) +
			  in(I,J-1,ek) + in(I-1,J+1,ek) ) * 0.25;

/// Twelve Edges
	out(I,0,0   ) = ( in(I,0,0) + in(I-1,0,0) ) * 0.5;
	out(I,0,ek+1) = ( in(I,0,ek) + in(I-1,0,ek) ) * 0.5;
	out(0   ,0,K) = ( in(0,0,K) + in(0  ,0,K-1) ) * 0.5;
	out(ei+1,0,K) = ( in(ei,0,K) + in(ei,0,K-1) ) * 0.5;
	out(I,ej+1,0   ) = ( in(I,ej,0) + in(I-1,ej,0) ) * 0.5;
	out(I,ej+1,ek+1) = ( in(I,ej,ek) + in(I-1,ej,ek) ) * 0.5;
	out(0,ej+1,K   ) = ( in(0,ej,K) + in(0,ej,K-1) ) * 0.5;
	out(ei+1,ej+1,K) = ( in(ei,ej,K) + in(ei,ej,K-1) ) * 0.5;
	out(0,J,0) = ( in(0,J,0) + in(0,J-1,0) ) * 0.5;
	out(0,J,ek+1) = ( in(0,J,ek) + in(0,J-1,ek) ) * 0.5;
	out(ei+1,J,0) = ( in(ei,J,0) + in(ei,J-1,0) ) * 0.5;
	out(ei+1,J,ek+1) = ( in(ei,J,ek) + in(ei,J-1,ek) ) * 0.5;
///
/// Eight corners
///
	out(0   ,0   ,0  ) = in(0,0,0);
	out(ei+1,0   ,0  ) = in(ei,0,0);
	out(0   ,ej+1,0  ) = in(0,ej,0);
	out(ei+1,ej+1,0  ) = in(ei,ej,0);
	out(0   ,0   ,ek+1) = in(0,0,ek);
	out(ei+1,0   ,ek+1) = in(ei,0,ek);
	out(0   ,ej+1,ek+1) = in(0,ej,ek);
	out(ei+1,ej+1,ek+1) = in(ei,ej,ek);	
    }    
}

/*!
 ***************************************************************************
 * Interpolate an scalarfield from center of volumes to nodes (1D,2D or 3D).
 * \verbatim
      in field (volumes)               out field (nodes)

5  x--x--+--x--+--x--+--x--x        o-----o-----o-----o-----o  4
   |     |     |     |     |        |     |     |     |     |
4  x  x  |  x  |  x  |  x  x        |     |     |     |     |
   |     |     |     |     |        |     |     |     |     |
   +-----+-----+-----+-----+        o-----o-----o-----o-----o  3   
   |     |     |     |     |        |     |     |     |     |
3  x  x  |  x  |  x  |  x  x        |     |     |     |     |
   |     |     |     |     |        |     |     |     |     |
   +-----+-----+-----+-----+  ===>  o-----o-----o-----o-----o  2   
   |     |     |     |     |        |     |     |     |     |
2  x  x  |  x  |  x  |  x  x        |     |     |     |     |
   |     |     |     |     |        |     |     |     |     |
   +-----+-----+-----+-----+        o-----o-----o-----o-----o  1   
   |     |     |     |     |        |     |     |     |     |
1  x  x  |  x  |  x  |  x  x        |     |     |     |     |
   |     |     |     |     |        |     |     |     |     |
0  x--x--+--x--+--x--+--x--x        o-----o-----o-----o-----o  0   

   0  1     2     3     4  5        0     1     2     3     4 

 input size : (vols_x, vols_y)     output size : (nodes_x, nodes_y)    
 Note: vols_? = nodes_? + 1, for ? = x, y, z \endverbatim
 ***************************************************************************
 *  \author Luis M. de la Cruz [ Sun Mar 23 15:33:45 GMT 2008 ]
 ***************************************************************************
 */
template<class Tfield>
void NumUtils::interpolateToNodes(Tfield &out, Tfield &in) 
{
    int bi = out.lbound(firstDim) + 1, ei = out.ubound(firstDim) - 1, 
	bj = out.lbound(secondDim)+ 1, ej = out.ubound(secondDim)- 1,
	bk = out.lbound(thirdDim) + 1, ek = out.ubound(thirdDim) - 1;
    
    Range I(bi,ei), J(bj,ej), K(bk,ek);

    if (out.rank() == 1) { 
	out(I) = ( in(I) + in(I+1) ) * 0.5;
	out(0) = in(0,0);
	out(ei+1) = in(ei+2);
    }
    if (out.rank() == 2) {
//
// Internal points    
//
	out(I,J) = ( in(I,J) + in(I+1,J) + in(I,J+1) + in(I+1,J+1) ) * 0.25;
//	
// Boundaries
//
	out(I   ,0   ) = ( in(I   ,0   ) + in(I+1 ,0   ) ) * 0.5;
	out(I   ,ej+1) = ( in(I   ,ej+2) + in(I+1 ,ej+2) ) * 0.5;
	out(0   ,J   ) = ( in(0   ,J   ) + in(0   ,J+1 ) ) * 0.5;
	out(ei+1,J   ) = ( in(ei+2,J   ) + in(ei+2,J+1 ) ) * 0.5;
//	
// Corners
//
	out(0   ,0   ) = in(0   ,0   );
	out(ei+1,0   ) = in(ei+2,0   );
	out(0   ,ej+1) = in(0   ,ej+2);
	out(ei+1,ej+1) = in(ei+2,ej+2);
    }
    if (out.rank() == 3) {
//
// Internal points    
//
	out(I,J,K) = ( in(I,J  ,K  ) + in(I+1,J  ,K  ) + 
		       in(I,J+1,K  ) + in(I+1,J+1,K  ) +
		       in(I,J  ,K+1) + in(I+1,J  ,K+1) + 
		       in(I,J+1,K+1) + in(I+1,J+1,K+1) ) * 0.125;
///	
/// Boundaries
///
/// Bot
	out(I,0   ,K) = ( in(I,0,K  ) + in(I+1,0,K  ) +
			  in(I,0,K+1) + in(I+1,0,K+1) ) * 0.25;
/// Top
	out(I,ej+1,K) = ( in(I,ej+2,K  ) + in(I+1,ej+2,K  ) +
			  in(I,ej+2,K+1) + in(I+1,ej+2,K+1) ) * 0.25;
/// Left
	out(0   ,J,K) = ( in(0,J,K  ) + in(0,J+1,K  ) +
			  in(0,J,K+1) + in(0,J+1,K+1) ) * 0.25;
/// Right
	out(ei+1,J,K) = ( in(ei+2,J,K  ) + in(ei+2,J+1,K  ) +
			  in(ei+2,J,K+1) + in(ei+2,J+1,K+1) ) * 0.25;
/// Back
	out(I,J,0   ) = ( in(I,J  ,0) + in(I+1,J  ,0) +
			  in(I,J+1,0) + in(I+1,J+1,0) ) * 0.25;
/// Front
	out(I,J,ek+1) = ( in(I,J  ,ek+2) + in(I+1,J  ,ek+2) +
			  in(I,J+1,ek+2) + in(I+1,J+1,ek+2) ) * 0.25;

/// Twelve Edges
	out(I,0,0   ) = ( in(I,0,0   ) + in(I+1,0,0   ) ) * 0.5;
	out(I,0,ek+1) = ( in(I,0,ek+2) + in(I+1,0,ek+2) ) * 0.5;
	out(0   ,0,K) = ( in(0   ,0,K) + in(0   ,0,K+1) ) * 0.5;
	out(ei+1,0,K) = ( in(ei+2,0,K) + in(ei+2,0,K+1) ) * 0.5;
	out(I,ej+1,0   ) = ( in(I,ej+2,0) + in(I+1,ej+2,0) ) * 0.5;
	out(I,ej+1,ek+1) = ( in(I,ej+2,ek+2) + in(I+1,ej+2,ek+2) ) * 0.5;
	out(0,ej+1,K   ) = ( in(0,ej+2,K) + in(0,ej+2,K+1) ) * 0.5;
	out(ei+1,ej+1,K) = ( in(ei+2,ej+2,K) + in(ei+2,ej+2,K+1) ) * 0.5;
	out(0,J,0) = ( in(0,J,0) + in(0,J+1,0) ) * 0.5;
	out(0,J,ek+1) = ( in(0,J,ek+2) + in(0,J+1,ek+2) ) * 0.5;
	out(ei+1,J,0) = ( in(ei+2,J,0) + in(ei+2,J+1,0) ) * 0.5;
	out(ei+1,J,ek+1) = ( in(ei+2,J,ek+2) + in(ei+2,J+1,ek+2) ) * 0.5;
///
/// Eight corners
///
	out(0   ,0   ,0  ) = in(0   ,0   ,0  );
	out(ei+1,0   ,0  ) = in(ei+2,0   ,0  );
	out(0   ,ej+1,0  ) = in(0   ,ej+2,0  );
	out(ei+1,ej+1,0  ) = in(ei+2,ej+2,0  );
	out(0   ,0   ,ek+1) = in(0   ,0   ,ek+2);
	out(ei+1,0   ,ek+1) = in(ei+2,0   ,ek+2);
	out(0   ,ej+1,ek+1) = in(0   ,ej+2,ek+2);
	out(ei+1,ej+1,ek+1) = in(ei+2,ej+2,ek+2);	
    }
    
}


/*!
 ***************************************************************************
 * Interpolate an scalarfield from center of volumes to nodes (1D,2D or 3D).
 * for Structured and non Uniform Meshes
 * \verbatim
      in field (volumes)               out field (nodes)

5  x--x--+--x--+--x--+--x--x        o-----o-----o-----o-----o  4
   |     |     |     |     |        |     |     |     |     |
4  x  x  |  x  |  x  |  x  x        |     |     |     |     |
   |     |     |     |     |        |     |     |     |     |
   +-----+-----+-----+-----+        o-----o-----o-----o-----o  3   
   |     |     |     |     |        |     |     |     |     |
3  x  x  |  x  |  x  |  x  x        |     |     |     |     |
   |     |     |     |     |        |     |     |     |     |
   +-----+-----+-----+-----+  ===>  o-----o-----o-----o-----o  2   
   |     |     |     |     |        |     |     |     |     |
2  x  x  |  x  |  x  |  x  x        |     |     |     |     |
   |     |     |     |     |        |     |     |     |     |
   +-----+-----+-----+-----+        o-----o-----o-----o-----o  1   
   |     |     |     |     |        |     |     |     |     |
1  x  x  |  x  |  x  |  x  x        |     |     |     |     |
   |     |     |     |     |        |     |     |     |     |
0  x--x--+--x--+--x--+--x--x        o-----o-----o-----o-----o  0   

   0  1     2     3     4  5        0     1     2     3     4 

 input size : (vols_x, vols_y)     output size : (nodes_x, nodes_y)    
 Note: vols_? = nodes_? + 1, for ? = x, y, z \endverbatim
 ***************************************************************************
 *  \author Luis M. de la Cruz [ Tue Jun 22 16:02:18 CDT 2010 ]
 ***************************************************************************
 * \todo 1D and 3D
 */
template<typename Tfield, typename Tmesh>
void NumUtils::interpolateToNodes(Tfield &out, Tfield &in, Tmesh &m) 
{
    int bi = out.lbound(firstDim) + 1, ei = out.ubound(firstDim) - 1, 
	bj = out.lbound(secondDim)+ 1, ej = out.ubound(secondDim)- 1,
	bk = out.lbound(thirdDim) + 1, ek = out.ubound(thirdDim) - 1;

    double alpha, beta;
    Range I(bi,ei), J(bj,ej), K(bk,ek);

    typedef typename TunaArray<double,1>::huge ScalarField1D;
    ScalarField1D dx = m.getDeltasVols(X), 
      dy = m.getDeltasVols(Y), dz = m.getDeltasVols(Y);

    if (out.rank() == 2) {
//
// Internal points
//
	for(int i = bi; i <= ei; ++i)
	    for(int j = bj; j <= ej; ++j) {
		alpha = dx(i) / ( dx(i) + dx(i+1) );
		beta  = dy(j) / ( dy(j) + dy(j+1) );	
		out(i,j) = ( in(i,j) * (1 - alpha) + 
			     in(i+1,j) * alpha ) * (1 - beta) +
		           ( in(i,j+1) * (1 - alpha) + 
			     in(i+1,j+1) * alpha ) * beta;
	    }
//
// Boundaries : Horizontal
//
	for(int i = bi; i <= ei; ++i){
	    alpha = dx(i) / ( dx(i) + dx(i+1) );
	    out(i,0) = in(i,0) * (1 - alpha) + in(i+1,0) * alpha;
	    out(i,ej+1) = in(i,ej+2) * (1 - alpha) + in(i+1,ej+2) * alpha;
	}
//
// Boundaries : Vertical
//
	for(int j = bj; j <= ej; ++j){
	    beta  = dy(j) / ( dy(j) + dy(j+1) );
	    out(0,j) = in(0,j) * (1 - beta) + in(0,j+1) * beta;
	    out(ei+1,j) = in(ei+2,j) * (1 - beta) + in(ei+2,j+1) * beta;
	}
//
//  Corners
//
	out(0   ,0   ) = in(0   ,0   );
	out(ei+1,0   ) = in(ei+2,0   );
	out(0   ,ej+1) = in(0   ,ej+2);
	out(ei+1,ej+1) = in(ei+2,ej+2);
    }
}


/*! 
 *  This function makes a shift of velocities from nodes
 *  to volume faces in x direction.
 */
template<class Tfield> 
void NumUtils::staggerFromNodesX(Tfield &output, Tfield &input)
{
    int bi = output.lbound(firstDim),    ei = output.ubound(firstDim);    
    int bj = output.lbound(secondDim)+1, ej = output.ubound(secondDim)-1;
    int bk = output.lbound(thirdDim)+1,  ek = output.ubound(thirdDim)-1;
    Range I(bi, ei), J(bj, ej), K(bk, ek);
    int Dim = output.rank();

    if (Dim == 1) {
      // In this case the size of input and output is the same.
      // The location of nodes and faces coincide in space.
        output(I) = input(I); 
    } else if (Dim == 2) {
	output(I,J) = ( input(I,J) + input(I,J-1) ) * 0.5;
	output(I,bj-1) = input(I,bj-1);
	output(I,ej+1) = input(I,ej); //! extravolumes are on boundary
    } else if (Dim == 3) {
	output(I,J,K) = ( input(I,J,K  ) + input(I,J-1,K  ) +
			  input(I,J,K-1) + input(I,J-1,K-1) ) * 0.25;
	output(I,bj-1,K   ) = input(I,bj-1,K   ); // bot wall
	output(I,ej+1,K   ) = input(I,ej  ,K   ); // top wall
	output(I,J   ,bk-1) = input(I,J   ,bk-1); // back wall
	output(I,J   ,ek+1) = input(I,J   ,ek  ); // front wall
	output(I,bj-1,bk-1) = input(I,bj-1,bk-1); // back-bot edge 
	output(I,ej+1,bk-1) = input(I,ej  ,bk-1); // back-top edge 
	output(I,ej+1,ek+1) = input(I,ej  ,ek  ); // front-top edge
	output(I,bj-1,ek+1) = input(I,bj-1,ek  ); // front-bot edge
    }
}

/*! 
 *  This function makes a shift of velocities from nodes
 *  to volume faces in y direction.
 */
template<class Tfield> 
void NumUtils::staggerFromNodesY(Tfield &output, Tfield &input)
{
    int bi = output.lbound(firstDim)+1, ei = output.ubound(firstDim)-1;
    int bj = output.lbound(secondDim), ej = output.ubound(secondDim);
    int bk = output.lbound(thirdDim)+1,  ek = output.ubound(thirdDim)-1;
    
    Range I(bi, ei), J(bj, ej), K(bk, ek);
    int Dim = output.rank();

    if (Dim == 2) {
	output(I,J) = ( input(I,J) + input(I-1,J) ) * 0.5;
	output(bi-1,J) = input(bi-1,J);
	output(ei+1,J) = input(ei,J); //! extravolumes are on boundary
    } else if (Dim == 3) {
	output(I,J,K) = ( input(I,J,K  ) + input(I-1,J,K  ) +
			  input(I,J,K-1) + input(I-1,J,K-1) ) * 0.25;
	output(bi-1,J,K   ) = input(bi-1,J,K   ); // left wall
	output(ei+1,J,K   ) = input(ei  ,J,K   ); // right wall
	output(I   ,J,bk-1) = input(I   ,J,bk-1); // back wall
	output(I   ,J,ek+1) = input(I   ,J,ek  ); // front wall
	output(bi-1,J,bk-1) = input(bi-1,J,bk-1); // back-left edge
	output(ei+1,J,bk-1) = input(ei  ,J,bk-1); // back-right edge
	output(bi-1,J,ek+1) = input(bi-1,J,ek  ); // front-left edge
	output(ei+1,J,ek+1) = input(ei  ,J,ek  ); // front-right edge
    }
}

/*! 
 *  This function makes a shift of velocities from nodes
 *  to volume faces in z direction.
 */
template<class Tfield> 
void NumUtils::staggerFromNodesZ(Tfield &output, Tfield &input)
{
    int bi = output.lbound(firstDim) +1, ei = output.ubound(firstDim)-1;
    int bj = output.lbound(secondDim)+1, ej = output.ubound(secondDim)-1;
    int bk = output.lbound(thirdDim), ek = output.ubound(thirdDim);    
    Range I(bi, ei), J(bj, ej), K(bk,ek);

    output(I,J,K) = ( input(I,J,K  ) + input(I-1,J,K) + 
		      input(I,J-1,K) + input(I-1,J-1,K) ) * 0.25;
    output(I   ,bj-1,K) = input(I   ,bj-1,K); // bot wall
    output(I   ,ej+1,K) = input(I   ,ej  ,K); // top wall
    output(bi-1,J   ,K) = input(bi-1,J   ,K); // left wall
    output(ei+1,J   ,K) = input(ei  ,J   ,K); // right wall
    output(bi-1,bj-1,K) = input(bi-1,bj-1,K); // left-bot edge
    output(bi-1,ej+1,K) = input(bi-1,ej  ,K); // left-top edge
    output(ei+1,bj-1,K) = input(ei  ,bj-1,K); // right-bot edge
    output(ei+1,ej+1,K) = input(ei  ,ej  ,K); // right-top edge
}

template<typename Tfield>
double NumUtils::calcErrorL1(const Tfield& phi1, const Tfield& phi2)
{
  Tfield ERROR (phi1 - phi2 );
  return NumUtils::normL1( ERROR ) / phi1.numElements();
}

template<typename Tfield>
double NumUtils::calcAbsoluteError_L1(const Tfield& phi1, const Tfield& phi2)
{
  Tfield ERROR (phi1 - phi2 );
  double norm = NumUtils::normL1( ERROR );
  //  cout << "\n RELATIVE = " << relative << "\t ABSOLUTE = " << norm << endl; 
  return norm;
}

template<typename Tfield>
double NumUtils::calcErrorL2(const Tfield& phi1, const Tfield& phi2)
{
  Tfield ERROR (phi1 - phi2 );
  double normL2 =  NumUtils::normL2( ERROR );
  //  return normL2 * normL2 / phi1.numElements(); // MSE
  return normL2 / sqrt(phi1.numElements()); // RMS
}

template<typename Tfield>
double NumUtils::calcAbsoluteError_L2(const Tfield& phi1, const Tfield& phi2)
{
  Tfield ERROR (phi1 - phi2 );
  double norm = NumUtils::normL2( ERROR );
  //  cout << "\n RELATIVE = " << relative << "\t ABSOLUTE = " << norm; 
  return norm;
}

template<typename Tfield>
double NumUtils::calcErrorLmax(const Tfield& phi1, const Tfield& phi2)
{
  Tfield ERROR (phi1 - phi2 );
  return NumUtils::normLmax( ERROR );
}

template<typename Tfield>
double NumUtils::calcAbsoluteError_Lmax(const Tfield& phi1, const Tfield& phi2)
{
  Tfield ERROR (phi1 - phi2);
  double norm = NumUtils::normLmax( ERROR );
  return norm;
}


template<typename Tfield>
double NumUtils::normL1(const Tfield &sf)
{
  double norm = 0;
  int bi, ei, bj, ej, bk, ek;
  int Dim = sf.rank();

  bi = sf.lbound(firstDim); ei = sf.ubound(firstDim);
  
  if (Dim == 1) {

    for(int i = bi; i <= ei; ++i) 
      norm += fabs ( sf(i) );
    
  } else if (Dim == 2) {
    
    bj = sf.lbound(secondDim); ej = sf.ubound(secondDim);    
    for(int i = bi; i <= ei; ++i)
      for(int j = bj; j <= ej; ++j)
	norm += fabs( sf(i,j) );
    
  } else if (Dim == 3) {
    
    bj = sf.lbound(secondDim); ej = sf.ubound(secondDim);
    bk = sf.lbound(thirdDim);  ek = sf.ubound(thirdDim);
  
    for(int k = bk; k <= ek; ++k)
      for(int i = bi; i <= ei; ++i)
	for(int j = bj; j <= ej; ++j)
	  norm += fabs( sf(i,j,k) );
  }
  return norm;
}

template<typename Tfield>
double NumUtils::normL2(const Tfield &sf) 
{
  double norm;
  int bi, ei, bj, ej, bk, ek;
  int Dim = sf.rank();

  bi = sf.lbound(firstDim); ei = sf.ubound(firstDim);

  norm = 0.0;
  if (Dim == 1) {

    for(int i = bi; i <= ei; ++i) 
      norm += ( sf(i) * sf(i) );

  } else if (Dim == 2) {

    bj = sf.lbound(secondDim); ej = sf.ubound(secondDim);    
    for(int i = bi; i <= ei; ++i)
      for(int j = bj; j <= ej; ++j)
	norm += ( sf(i,j) * sf(i,j) );

  } else if (Dim == 3) {

    bj = sf.lbound(secondDim); ej = sf.ubound(secondDim);
    bk = sf.lbound(thirdDim);  ek = sf.ubound(thirdDim);    
    for(int k = bk; k <= ek; ++k)
      for(int i = bi; i <= ei; ++i)
	for(int j = bj; j <= ej; ++j)
	  norm += ( sf(i,j,k) * sf(i,j,k) );
  }
  return sqrt (norm);
}


template<typename Tfield>
double NumUtils::normLmax(const Tfield &sf)
{
  double max = 0.0;
  int bi, ei, bj, ej, bk, ek;
  int Dim = sf.rank();

  bi = sf.lbound(firstDim); ei = sf.ubound(firstDim);

  if (Dim == 1) {

    for(int i = bi; i <= ei; ++i)
      if ( fabs( sf(i) ) > max ) max = fabs( sf(i) ); 

  } else if (Dim == 2) {

    bj = sf.lbound(secondDim); ej = sf.ubound(secondDim);    
    for(int i = bi; i <= ei; ++i)
      for(int j = bj; j <= ej; ++j)	
	if ( fabs( sf(i,j) ) > max ) max = fabs( sf(i,j) );

  } else if (Dim == 3) {

    bj = sf.lbound(secondDim); ej = sf.ubound(secondDim);
    bk = sf.lbound(thirdDim);  ek = sf.ubound(thirdDim);    
    for(int k = bk; k <= ek; ++k)
      for(int i = bi; i <= ei; ++i)
	for(int j = bj; j <= ej; ++j)
	  if ( fabs( sf(i,j,k) ) > max ) max = fabs ( sf(i,j,k) );
  }
  return max;
}

/*!
 * Interpolation from staggered mesh to the nodes of u-velocity.
 */
template<class Tfield>
void NumUtils::interpolateToNodesUs(Tfield &out, Tfield &in) 
{
    int bi = out.lbound(firstDim), ei = out.ubound(firstDim), 
	bj = out.lbound(secondDim)+ 1, ej = out.ubound(secondDim)- 1,
	bk = out.lbound(thirdDim) + 1, ek = out.ubound(thirdDim) - 1;
    
    Range I(bi,ei), J(bj,ej), K(bk,ek);

    if (out.rank() == 2) {
//
// Internal points
//
	out(I,J) = ( in(I,J) + in(I,J+1) ) * 0.5;	
//
// Boundaries 
// 
// Top
//
	out(I ,0 ) = in(I ,0 );
//
// Bottom
//
	out(I ,ej+1) = in(I ,ej+2);
    }
    if (out.rank() == 3) {
//
// Internal points
//
	out(I,J,K) = ( in(I,J,K  ) + in(I,J+1,K  ) + 
		       in(I,J,K+1) + in(I,J+1,K+1) ) * 0.25;
//
// Boundaries
//
// Bot
	out(I,0   ,K) = ( in(I,0   ,K) + in(I,0   ,K+1) ) * 0.5;	
// Top
	out(I,ej+1,K) = ( in(I,ej+2,K) + in(I,ej+2,K+1) ) * 0.5;
// Front
	out(I,J,0   ) = ( in(I,J,0   ) + in(I,J+1,0   ) ) * 0.5;
// Back
	out(I,J,ek+1) = ( in(I,J,ek+2) + in(I,J+1,ek+2) ) * 0.5;
    }    
}

/*!
 * Interpolation from staggered mesh to the nodes of v-velocity.
 */
template<class Tfield>
void NumUtils::interpolateToNodesVs(Tfield &out, Tfield &in) 
{
    int bi = out.lbound(firstDim) + 1, ei = out.ubound(firstDim) - 1, 
	bj = out.lbound(secondDim), ej = out.ubound(secondDim),
	bk = out.lbound(thirdDim) + 1, ek = out.ubound(thirdDim) - 1;
    
    Range I(bi,ei), J(bj,ej), K(bk,ek);

    if (out.rank() == 2) {
//
// Internal points
//
	out(I,J) = ( in(I,J) + in(I+1,J) ) * 0.5;	
//
// Boundaries (just right and left walls, and no corners)
//
	out(0   , J) = in(0   , J);
	out(ei+1, J) = in(ei+2, J);
    }
    if (out.rank() == 3) {
//
// Internal points
//
	out(I,J,K) = ( in(I,J,K  ) + in(I+1,J,K  ) + 
		       in(I,J,K+1) + in(I+1,J,K+1) ) * 0.25;
//
// Boundaries
//
// Left
	out(0   ,J,K) = ( in(0   ,J,K) + in(0   ,J,K+1) ) * 0.5;	
// Right
	out(ei+1,J,K) = ( in(ei+2,J,K) + in(ei+2,J,K+1) ) * 0.5;
// Front
	out(I,J,0   ) = ( in(I,J,0   ) + in(I+1,J,0   ) ) * 0.5;
// Back
	out(I,J,ek+1) = ( in(I,J,ek+2) + in(I+1,J,ek+2) ) * 0.5;
    }    
}

/*!
 * Interpolation from staggered mesh to the nodes of w-velocity.
 */
template<class Tfield>
void NumUtils::interpolateToNodesWs(Tfield &out, Tfield &in) 
{
    int bi = out.lbound(firstDim) + 1, ei = out.ubound(firstDim) - 1, 
	bj = out.lbound(secondDim)+ 1, ej = out.ubound(secondDim)- 1,
	bk = out.lbound(thirdDim), ek = out.ubound(thirdDim);
    
    Range I(bi,ei), J(bj,ej), K(bk,ek);

//
// Internal points
//
    out(I,J,K) = ( in(I,J  ,K) + in(I+1,J  ,K) + 
		   in(I,J+1,K) + in(I+1,J+1,K) ) * 0.25;
//
// Boundaries
//
// Left
    out(0   ,J,K) = ( in(0   ,J,K) + in(0   ,J+1,K) ) * 0.5;	
// Right
    out(ei+1,J,K) = ( in(ei+2,J,K) + in(ei+2,J+1,K) ) * 0.5;
// Bot
    out(I,0   ,K) = ( in(I,0   ,K) + in(I+1,0   ,K) ) * 0.5;
// Top
    out(I,ej+1,K) = ( in(I,ej+2,K) + in(I+1,ej+2,K) ) * 0.5;    
}


template<class Tfield> 
bool NumUtils::interX(Tfield &output, Tfield &input)
//!----------------------------------------------------------------------
//!        output field                       input field 
//!
//!  4  >   >   >   >   >                4  >   >   >   >
//!       +---+---+---+                     +---+---+---+    
//!  3  > | > | > | > | >                3  > * > * > * >
//!       +---+---+---+                     +---+---+---+        
//!  2  > | > | > | > | >     interX     2  > * > * > * >
//!       +---+---+---+     <=========      +---+---+---+    
//!  1  > | > | > | > | >                1  > * > * > * >
//!       +---+---+---+                     +---+---+---+    
//!  0  >   >   >   >   >                0  >   >   >   >
//!                                   
//!     0   1   2   3   4                   0   1   2   3  
//!
//!        Centred                            Staggered
//!        -------                            ---------
//! Mesh Nodes(+): nx = 3; ny = 3;       Mesh Nodes(+): nx = 3; ny = 3;        
//! Vols. uvel(>): vx = 4; vy = 4;       Vols. uvel(>): vsx = 3; vsy = 4;
//!
//! output size : (vx, vy)                input size : (vsx, vsy)
//!              
//!----------------------------------------------------------------------
{
    int bi = output.lbound(firstDim)+1,  ei = output.ubound(firstDim)-1,
	bj = output.lbound(secondDim),  ej = output.ubound(secondDim);
    Range I(bi, ei), J(bj, ej);
    int Dim = output.rank();
    
    if (Dim == 2) {
	output(I,J) = ( input(I-1,J) + input(I,J) ) * 0.5;
	output(bi-1, J) = input(bi-1,J);
	output(ei+1, J) = input(ei+1,J); //! extravolumes are on boundary
    } else if (Dim == 3) {
	int bk = output.lbound(thirdDim),  ek = output.ubound(thirdDim);
	Range K(bk, ek);
	output(I,J,K) = ( input(I-1,J,K) + input(I,J,K) ) * 0.5;
	output(bi-1,J,K) = input(bi-1,J,K);
	output(ei+1,J,K) = input(ei+1,J,K); //! extravolumes are on boundary
    }
    
    return 1;
}

template<class Tfield> 
bool NumUtils::interY(Tfield &output, Tfield &input)
//!----------------------------------------------------------------------
//!        output field                        input field 
//!
//!  4  ^   ^   ^   ^   ^                   
//!       +---+---+---+                  3 ^ +-^-+-^-+-^-+ ^  
//!  3  ^ | ^ | ^ | ^ | ^                    | * | * | * |
//!       +---+---+---+                  2 ^ +-^-+-^-+-^-+ ^      
//!  2  ^ | ^ | ^ | ^ | ^      interX        | * | * | * |
//!       +---+---+---+     =========>   1 ^ +-^-+-^-+-^-+ ^  
//!  1  ^ | ^ | ^ | ^ | ^                    | * | * | * |
//!       +---+---+---+                  0 ^ +-^-+-^-+-^-+ ^  
//!  0  ^   ^   ^   ^   ^                  0   1   2   3   4 
//!     0   1   2   3   4                  
//!
//!        Centred                            Staggered
//!        -------                            ---------
//! Mesh Nodes(+): nx = 3; ny = 3;       Mesh Nodes(+): nx = 3; ny = 3;        
//! Vols. uvel(^): vx = 4; vy = 4;       Vols. uvel(^): vsx = 4; vsy = 3;
//!
//! output size : (vx, vy)               input size : (vsx, vsy)
//!              
//!----------------------------------------------------------------------
{
    int bi = output.lbound(firstDim),  ei = output.ubound(firstDim),
	bj = output.lbound(secondDim)+1, ej = output.ubound(secondDim)-1;
    Range I(bi, ei), J(bj, ej);
    int Dim = output.rank();

    if (Dim == 2) {
	output(I,J) = ( input(I,J-1) + input(I,J) ) * 0.5;
	output(I,bj-1) = input(I,bj-1);
	output(I,ej+1) = input(I,ej+1); //! extravolumes are on boundary
    } else if (Dim == 3) {
	int bk = output.lbound(thirdDim),  ek = output.ubound(thirdDim);
	Range K(bk, ek);
	output(I,J,K) = ( input(I,J-1,K) + input(I,J,K) ) * 0.5;
	output(I,bj-1,K) = input(I,bj-1,K);
	output(I,ej+1,K) = input(I,ej+1,K); //! extravolumes are on boundary
    }
    
    return 1;
}


template<class Tfield> 
bool NumUtils::interZ(Tfield &output, Tfield &input)
{
    int bi = output.lbound(firstDim),  ei = output.ubound(firstDim),
	bj = output.lbound(secondDim), ej = output.ubound(secondDim),
	bk = output.lbound(thirdDim)+1,  ek = output.ubound(thirdDim)-1;
    Range I(bi, ei), J(bj, ej), K(bk, ek);
    int Dim = output.rank();
    
    output(I,J,K) = ( input(I,J,K-1) + input(I,J,K) ) * 0.5;
    output(I,J,bk-1) = input(I,J,bk-1);
    output(I,J,ek+1) = input(I,J,ek+1); //! extravolumes are on boundary

    return 1;
}

/*!
 * Calculate de Nusselt number.
 * \param T the temperature scalar field
 * \param dx delta in x-axis
 * \param dy delta in y-axis
 * \param mon i-column of the mesh
 * \par Description
 *      This function calculate the Nusselt number in a vertical line 
 *      of the 2D rectangular domain where the temperature was calculated.
 *      The integer mon defines the i-column where the nusselt will be 
 *      calculated. For example, if this value is zero, then the nusselt
 *      will be calculated in the left wall of the rectangular domain.
 *      The Nusselt number is defined as:
 * \[
 * \]
 */
template<class Tfield>
double NumUtils::calcNusseltVert(Tfield &T, double dx, double dy, int mon) 
{
    int bj = T.lbound(secondDim), ej = T.ubound(secondDim);
    double dy_dx = dy / dx, nu = 0;
    
    for(int j = bj; j <= ej; ++j)
      nu -= ( T(mon+1, j) - T(mon, j) ) * dy_dx;
    return nu;
}

/*!
 * The same as the above function but for horizontal lines.
 */
template<class Tfield>
double NumUtils::calcNusseltHoriz(Tfield &T, double dx, double dy, int mon)
{
    int bi = T.lbound(firstDim), ei = T.ubound(firstDim)-1;
    double dx_dy = dx / dy, nu = 0;
    
    for(int i = bi; i <= ei; ++i)
      nu -= ( T(i, mon+1) - T(i, mon) ) * dx_dy;
    
    return nu;
}


/*!
 * Calculate de Nusselt number.
 * \param T the temperature scalar field
 * \param dx delta in x-axis
 * \param dy delta in y-axis
 * \param dz delta in z-axis
 * \param mon i-column of the mesh
 * \par Description
 *      This function calculate the Nusselt number in a vertical plane,
 *      parallel to yz plane, of the 3D rectangular domain where the 
 *      temperature was calculated.
 *      The integer mon defines the i-column where the nusselt will be 
 *      calculated. For example, if this value is zero, then the nusselt
 *      will be calculated in the left wall of the rectangular domain.
 *      The Nusselt number is defined as:
 * \[ 
 * \]
 */
template<class Tfield>
double NumUtils::calcNusseltVert(Tfield &T, 
				 double dx, double dy, double dz, int mon) 
{
    int bj = T.lbound(secondDim), ej = T.ubound(secondDim);
    int bk = T.lbound(thirdDim), ek = T.ubound(thirdDim);

    double dyz_dx = dy * dz / dx, nuc = 0;

    for(int j = bj; j <= ej; ++j) 
	for(int k = bk; k <= ek; ++k) 
	    nuc -= ( T(mon+1,j,k) - T(mon,j,k) ) * dyz_dx;

    return nuc;
}

template<class Tfield>
double NumUtils::calcNusseltHoriz(Tfield &T, 
				  double dx, double dy, double dz, int mon) 
{
    int bi = T.lbound(firstDim), ei = T.ubound(firstDim);
    int bk = T.lbound(thirdDim), ek = T.ubound(thirdDim);

    double dxz_dy = dx * dz / dy, nuc = 0;

    for(int i = bi; i <= ei; ++i) 
	for(int k = bk; k <= ek; ++k) 
	    nuc -= ( T(i,mon+1,k) - T(i,mon,k) ) * dxz_dy;

    return nuc;
}

/***----------------------------------------------------------------------***

/// This stagger function makes a shift of velocities from volume
/// centers to volume faces.
template<class Tfield> 
bool NumUtils::staggerX(Tfield &output, Tfield &input)
//!----------------------------------------------------------------------
//!        input field                       output field 
//!
//!  4  >   >   >   >   >                4  >   >   >   >
//!       +---+---+---+                     +---+---+---+    
//!  3  > | > | > | > | >                3  > * > * > * >
//!       +---+---+---+                     +---+---+---+        
//!  2  > | > | > | > | >    staggerX    2  > * > * > * >
//!       +---+---+---+     =========>      +---+---+---+    
//!  1  > | > | > | > | >                1  > * > * > * >
//!       +---+---+---+                     +---+---+---+    
//!  0  >   >   >   >   >                0  >   >   >   >
//!                                   
//!     0   1   2   3   4                   0   1   2   3  
//!
//!        Centred                            Staggered
//!        -------                            ---------
//! Mesh Nodes(+): nx = 3; ny = 3;       Mesh Nodes(+): nx = 3; ny = 3;        
//! Vols. uvel(>): vx = 4; vy = 4;       Vols. uvel(>): vsx = 3; vsy = 4;
//!
//! input size : (vx, vy)                output size : (vsx, vsy)
//!              
//!----------------------------------------------------------------------
{
    int bi = output.lbound(firstDim)+1, ei = output.ubound(firstDim)-1,
	bj = output.lbound(secondDim), ej = output.ubound(secondDim);

    Range I(bi, ei), J(bj, ej);
    int Dim = output.rank();
    
    if (Dim == 2) {
	output(I,J) = ( input(I,J) + input(I+1,J) ) * 0.5;
	output(bi-1, J) = input(bi-1,J);
	output(ei+1, J) = input(ei+1,J); //! extravolumes are on boundary
    } else if (Dim == 3) {
	int bk = output.lbound(thirdDim),  ek = output.ubound(thirdDim);
	Range K(bk, ek);
	output(I,J,K) = ( input(I,J,K) + input(I+1,J,K) ) * 0.5;
	output(bi-1,J,K) = input(bi-1,J,K);
	output(ei+1,J,K) = input(ei+1,J,K); //! extravolumes are on boundary
    }
    
    return 1;
}

/// This stagger function makes a shift of velocities from nodes
/// to volume faces, for non-uniform meshes.
template<class Tfield, class Tmesh> 
bool NumUtils::staggerX(Tmesh &m, Tfield &output, Tfield &input)
{
    int bi = output.lbound(firstDim), ei = output.ubound(firstDim),
	bj = output.lbound(secondDim)+1, ej = output.ubound(secondDim)-1;
    Range I(bi, ei), J(bj, ej);
    int Dim = output.rank();

    double beta, gamma;
    ScalarField1D dy = m.getDeltaY_ns(), dz = m.getDeltaZ_fb();

    if (Dim == 2) {
	for(int i = bi; i <= ei; ++i)
	    for(int j = bj; j <= ej; ++j) {
		beta = dy(j) / ( dy(j) + dy(j-1) );
		output(i,j) = input(i,j) * (1 - beta) + input(i,j-1) * beta;
	    }
	output(I,bj-1) = input(I,bj-1);
	output(I,ej+1) = input(I,ej); 
    } else if (Dim == 3) {
	int bk = output.lbound(thirdDim)+1,  ek = output.ubound(thirdDim)-1;
	Range K(bk, ek);
	for(int i = bi; i <= ei; ++i)
	    for(int j = bj; j <= ej; ++j)
		for(int k = bk; k <= ek; ++k){
		    beta = dy(j) / ( dy(j) + dy(j-1) );
		    gamma = dz(k) / ( dz(k) + dz(k-1) );
		    output(i,j) = ( input(i,j  ,k) * (1 - beta) + 
				    input(i,j-1,k) * beta ) * (1 - gamma) +
                                  ( input(i,j  ,k-1) * (1 - beta) + 
				    input(i,j-1,k-1) * beta ) * gamma;
		}
	output(I,bj-1,K   ) = input(I,bj-1,K   ); // bot wall
	output(I,ej+1,K   ) = input(I,ej  ,K   ); // top wall
	output(I,J   ,bk-1) = input(I,J   ,bk-1); // back wall
	output(I,J   ,ek+1) = input(I,J   ,ek  ); // front wall
	output(I,bj-1,bk-1) = input(I,bj-1,bk-1); // back-bot edge
	output(I,ej+1,bk+1) = input(I,ej  ,bk  ); // back-top edge 
	output(I,ej+1,ek+1) = input(I,ej  ,ek  ); // front-top edge
	output(I,bj-1,ek+1) = input(I,bj-1,ek  ); // front-bot edge
    }

    return 1;
}


template<class Tfield> 
bool NumUtils::staggerY(Tfield &output, Tfield &input)
//!----------------------------------------------------------------------
//!        input field                       output field 
//!
//!  4  ^   ^   ^   ^   ^                   
//!       +---+---+---+                  3 ^ +-^-+-^-+-^-+ ^  
//!  3  ^ | ^ | ^ | ^ | ^                    | * | * | * |
//!       +---+---+---+                  2 ^ +-^-+-^-+-^-+ ^      
//!  2  ^ | ^ | ^ | ^ | ^    staggerX        | * | * | * |
//!       +---+---+---+     =========>   1 ^ +-^-+-^-+-^-+ ^  
//!  1  ^ | ^ | ^ | ^ | ^                    | * | * | * |
//!       +---+---+---+                  0 ^ +-^-+-^-+-^-+ ^  
//!  0  ^   ^   ^   ^   ^                  0   1   2   3   4 
//!     0   1   2   3   4                  
//!
//!        Centred                            Staggered
//!        -------                            ---------
//! Mesh Nodes(+): nx = 3; ny = 3;       Mesh Nodes(+): nx = 3; ny = 3;        
//! Vols. uvel(^): vx = 4; vy = 4;       Vols. uvel(^): vsx = 4; vsy = 3;
//!
//! input size : (vx, vy)                output size : (vsx, vsy)
//!              
//!----------------------------------------------------------------------
{
    int bi = output.lbound(firstDim),  ei = output.ubound(firstDim),
	bj = output.lbound(secondDim)+1, ej = output.ubound(secondDim)-1;
    Range I(bi, ei), J(bj, ej);
    int Dim = output.rank();

    if (Dim == 2) {
	output(I,J) = ( input(I,J) + input(I,J+1) ) * 0.5;
	output(I,bj-1) = input(I,bj-1);
	output(I,ej+1) = input(I,ej+1); //! extravolumes are on boundary
    } else if (Dim == 3) {
	int bk = output.lbound(thirdDim),  ek = output.ubound(thirdDim);
	Range K(bk, ek);
	output(I,J,K) = ( input(I,J,K) + input(I,J+1,K) ) * 0.5;
	output(I,bj-1,K) = input(I,bj-1,K);
	output(I,ej+1,K) = input(I,ej+1,K); //! extravolumes are on boundary
    }
    
    return 1;
}

template<class Tfield, class Tmesh> 
bool NumUtils::staggerY(Tmesh &m, Tfield &output, Tfield &input)
{
    int bi = output.lbound(firstDim)+1,  ei = output.ubound(firstDim)-1,
	bj = output.lbound(secondDim), ej = output.ubound(secondDim);
    Range I(bi, ei), J(bj, ej);
    int Dim = output.rank();
    double alpha, gamma;
    
    ScalarField1D dx = m.getDeltaX_ew(), dz = m.getDeltaZ_fb();

    if (Dim == 2) {
	for(int i = bi; i <= ei; ++i)
	    for(int j = bj; j <= ej; ++j) {
		alpha = dx(i) / ( dx(i) + dx(i-1) );
		output(i,j) = input(i,j) * (1 - alpha) + input(i-1,j) * alpha;
	    }
	output(bi-1,J) = input(bi-1,J);
	output(ei+1,J) = input(ei,J); 

    } else if (Dim == 3) {
	int bk = output.lbound(thirdDim)+1,  ek = output.ubound(thirdDim)-1;
	Range K(bk, ek);

	for(int i = bi; i <= ei; ++i)
	    for(int j = bj; j <= ej; ++j)
		for(int k = bk; k <= ek; ++k) {
		    alpha = dx(i) / ( dx(i) + dx(i-1) );
		    gamma = dz(k) / ( dz(k) + dz(k-1) );
		    output(i,j,k) = ( input(i  ,j,k) * (1 - alpha) + 
				      input(i-1,j,k) * alpha ) * (1 - gamma) +
			            ( input(i  ,j,k-1) * (1 - alpha) +
			              input(i-1,j,k-1) * alpha ) * gamma;
		}
	output(bi-1,J,K   ) = input(bi-1,J,K   ); // left wall
	output(ei+1,J,K   ) = input(ei  ,J,K   ); // right wall
	output(I   ,J,bk-1) = input(I   ,J,bk-1); // back wall
	output(I   ,J,ek+1) = input(I   ,J,ek  ); // front wall
	output(bi-1,J,bk-1) = input(bi-1,J,bk-1); // back-left edge
	output(ei+1,J,bk-1) = input(ei  ,J,bk-1); // back-right edge
	output(bi-1,J,ek+1) = input(bi-1,J,ek  ); // front-left edge
	output(ei+1,J,ek+1) = input(ei  ,J,ek  ); // front-right edge
    }
    
    return 1;
}


template<class Tfield> 
bool NumUtils::staggerZ(Tfield &output, Tfield &input)
{
    int bi = output.lbound(firstDim),  ei = output.ubound(firstDim),
	bj = output.lbound(secondDim), ej = output.ubound(secondDim),
	bk = output.lbound(thirdDim)+1,  ek = output.ubound(thirdDim)-1;
    Range I(bi, ei), J(bj, ej), K(bk, ek);
    int Dim = output.rank();

    output(I,J,K) = ( input(I,J,K) + input(I,J,K+1) ) * 0.5;
    output(I,J,bk-1) = input(I,J,bk-1);
    output(I,J,ek+1) = input(I,J,ek+1); //! extravolumes are on boundary

    return 1;
}

template<class Tfield, class Tmesh> 
bool NumUtils::staggerZ(Tmesh &m, Tfield &output, Tfield &input)
{
    int bi = output.lbound(firstDim),  ei = output.ubound(firstDim),
	bj = output.lbound(secondDim), ej = output.ubound(secondDim),
	bk = output.lbound(thirdDim)+1,  ek = output.ubound(thirdDim)-1;
    Range I(bi, ei), J(bj, ej), K(bk, ek);
    int Dim = output.rank();
    double alpha, beta;
    
    ScalarField1D dy = m.getDeltaY_fb(), dz = m.getDeltaZ_fb();

    for(int i = bi; i <= ei; ++i)
	for(int j = bj; j <= ej; ++j)
	    for(int k = bk; k <= ek; ++k) {
		alpha = dx(i) / ( dx(i) + dx(i-1) );
		beta  = dy(j) / ( dy(j) + dy(j-1) );
		output(i,j,k) = ( input(i  ,j  ,k) * (1 - alpha) + 
				  input(i-1,j  ,k) * alpha ) * (1 - beta) +
		                ( input(i  ,j-1,k) * (1 - alpha) +
                                  input(i-1,j-1,k) * alpha ) * beta;
	    }
    output(I   ,bj-1,K) = input(I   ,bj-1,K); // bot wall
    output(I   ,ej+1,K) = input(I   ,ej  ,K); // top wall
    output(bi-1,J   ,K) = input(bi-1,J   ,K); // left wall
    output(ei+1,J   ,K) = input(ei  ,J   ,K); // right wall
    output(bi-1,bj-1,K) = input(bi-1,bj-1,K); // left-bot edge
    output(bi-1,ej+1,K) = input(bi-1,ej  ,K); // left-top edge
    output(ei+1,bj-1,K) = input(ei  ,bj-1,K); // right-bot edge
    output(ei+1,ej+1,K) = input(ei  ,ej  ,K); // right-top edge

    return 1;
}



template<class Tfield>
double NumUtils::getSorsum(Tfield &u, Tfield &v, double dx, double dy) { 
    double sorsum = 0.0;

    int bi = 1, bj = 1;
    int ei = u.ubound(firstDim);
    int ej = u.ubound(secondDim) - 1;

    for(int i = bi; i <= ei; ++i ) 
	for(int j = bj; j <= ej; ++j)
	    sorsum += ( - u(i+1,j) * dy + u(i,j) * dy 
			- v(i,j+1) * dx + v(i,j) * dx );
    return sorsum; 
}

template<class Tfield>
double NumUtils::getSorsum(Tfield &u, Tfield &v, Tfield &w, 
		 double dx, double dy, double dz) { 
    double dyz = dy * dz, dxz = dx * dz, dxy = dx * dy;
    double sorsum = 0.0;

    int bi = 1, bj = 1, bk = 1;
    int ei = u.ubound(firstDim);
    int ej = u.ubound(secondDim) - 1;
    int ek = u.ubound(thirdDim) - 1;

    for(int i = bi; i <= ei; ++i )
	for(int j = bj; j <= ej; ++j)
	    for(int k = bk; k <= ek; ++k)
		sorsum += ( - u(i+1,j,k) * dyz + u(i,j,k) * dyz 
			    - v(i,j+1,k) * dxz + v(i,j,k) * dxz
			    - w(i,j,k+1) * dxy + w(i,j,k) * dxy );
    return sorsum; 
}


template<class Tfield>
double NumUtils::getResorm(Tfield &u, Tfield &v, double dx, double dy) { 
    double resorm = 0.0, smp = 0.0;
 
    int bi = 1, bj = 1;
    int ei = u.ubound(firstDim);
    int ej = u.ubound(secondDim) - 1;

    for(int i = bi; i <= ei; ++i ) {
	for(int j = bj; j <= ej; ++j)
	    smp = ( - u(i+1,j) * dy + u(i,j) * dy 
		    - v(i,j+1) * dx + v(i,j) * dx );
	resorm += fabs(smp);
    }

    return resorm;
}

template<class Tfield>
double NumUtils::getResorm(Tfield &u, Tfield &v, Tfield &w, 
		 double dx, double dy, double dz) { 
    double dyz = dy * dz, dxz = dx * dz, dxy = dx * dy;
    double resorm = 0.0, smp = 0.0;
 
    int bi = 1, bj = 1, bk = 1;
    int ei = u.ubound(firstDim);
    int ej = u.ubound(secondDim) - 1;
    int ek = u.ubound(thirdDim) - 1;   
    

    for(int i = bi; i <= ei; ++i ) 
	for(int j = bj; j <= ej; ++j) 
	    for(int k = bk; k <= ek; ++k)
	    {
		smp = ( - u(i+1,j,k) * dyz + u(i,j,k) * dyz 
			- v(i,j+1,k) * dxz + v(i,j,k) * dxz
			- w(i,j,k+1) * dxy + w(i,j,k) * dxy );
		resorm += fabs(smp);
	    }
    return resorm; 
}
***----------------------------------------------------------------------***/



/***----------------------------------------------------------------------***
template<class Tfield>
double NumUtils::calcNusselt(Tfield &T, double dx, double dy) 
{
    int bi = T.lbound(firstDim), ei = T.ubound(firstDim)-1;
    int bj = T.lbound(secondDim), ej = T.ubound(secondDim);
    double dy_dx = dy / dx, nuc = 0;
    int i;
    
    for(i = bi; i <= ei; ++i)
	for(int j = bj; j <= ej; ++j) 
	    nuc -= ( T(i+1,j) - T(i,j) ) * dy_dx;
    
    return nuc/i;
}
***----------------------------------------------------------------------***/


  /*
template<class Tfield, class Tmesh>
void NumUtils::interpolateToNodesUs(Tmesh &m, Tfield &out, Tfield &in) 
{
    int bi = out.lbound(firstDim), ei = out.ubound(firstDim),
	bj = out.lbound(secondDim)+ 1, ej = out.ubound(secondDim)- 1,
	bk = out.lbound(thirdDim) + 1, ek = out.ubound(thirdDim) - 1;
    double alpha, beta;
    ScalarField1D dy = m.getDeltaY();
    
    Range I(bi,ei), J(bj,ej), K(bk,ek);

    if (out.rank() == 2) {
//
// Internal points
//
	for(int i = bi; i <= ei; ++i)
	    for(int j = bj; j <= ej; ++j) {
		beta  = dy(j) / ( dy(j) + dy(j+1) );		
		out(i,j) = in(i,j) * (1.0 - beta) + in(i,j+1) * beta;
	    }
//
// Boundaries
//
	out(I ,0   ) = in(I ,0   );
	out(I ,ej+1) = in(I ,ej+2);
    }
    
}
  


template<class Tfield, class Tmesh>
void NumUtils::interpolateToNodesVs(Tmesh &m, Tfield &out, Tfield &in) 
{
    int bi = out.lbound(firstDim) + 1, ei = out.ubound(firstDim) - 1,
	bj = out.lbound(secondDim), ej = out.ubound(secondDim),
	bk = out.lbound(thirdDim) + 1, ek = out.ubound(thirdDim) - 1;
    double alpha, beta;
    ScalarField1D dx = m.getDeltaX();
    
    Range I(bi,ei), J(bj,ej), K(bk,ek);

    if (out.rank() == 2) {
//
// Internal points
//
	for(int i = bi; i <= ei; ++i)
	    for(int j = bj; j <= ej; ++j) {
		alpha = dx(i) / ( dx(i) + dx(i+1) );
		out(i,j) = in(i,j) * (1 - alpha) + in(i+1,j) * alpha;
	    }
//
// Boundaries
//	
	out(0   ,J) = in(0   ,J);
	out(ei+1,J) = in(ei+2,J);
    }
    
}


*/


} // Tuna namespace

